espHb = c()
espGb = c()
i = 1
for (b in seq(0,2,.01)){
espHb[i] = integrate(
function(x){
#Valores xi * f(x)
return(((x - b)^2)*dexp(x,2));
},
-Inf,
Inf)$value
espGb[i] = integrate(
function(x){
#Valores xi * f(x)
return((abs(x - b))*dexp(x,2));
},
-Inf,
Inf)$value
i = i+1
}
plot(seq(0,2,.01),
espHb,
lwd=3,
type='l',
col=2,
xlab='b',
ylab='Esperanza',
main='Esperanza en función de b',
)
lines(seq(0,2,.01),
espGb,
lwd=3,
type='l',
col=3
)
legend(x = "topleft",
legend = c("E((X-b)^2)", "E(|X-b|)"),
lty = c(1, 1),
col = c(2, 3)
)
axis(1,
at = seq(0,2,.1))
Consigna:
Determinar aproximadamente los valores b que minimizan la funcion H(b) y G(b) respectivamente
Respuesta:
Los valores que minimizan la funciones
parecieran rondar los:
0.5 para la función E((X-b)^2)
0.3 para la función E(|X-b|)
Consigna:
Tomar una muestra aleatoria de 10000 v.a. independientes exponenciales con λ = 2 y calcular i. la media(o promedio). ii. la mediana
va = rexp(10000, 2)
paste('Media: ', round(mean(va),4)); paste('Mediana: ', round(median(va),4))
## [1] "Media: 0.5008"
## [1] "Mediana: 0.3472"
Consigna:
Comparar la media anterior con b∗H, y la mediana con b∗G. ¿Son parecidos? Interprete.
Respuesta:
Los valores hallados son extremadamente
similares a los valores de b que minimizan las funciones del apartado
“a”. Más específicamente:
la media minimiza la función Hb, donde las
distancias se elevan al cuadrado.
Esto sucede porque, en particular,
elevar las distancias al cuadrado genera el efecto de una penalización
que crece a medida que las distancias crecen.
Restando la media
logramos obtener el vector que tiene las distancias “extremas” más
pequeñas, evitando en cuanto se pueda valores extremos.
la mediana minimiza la función Gb, donde se
calcula el módulo de las distancias.
Esto sucede porque, al no
elevar al cuadrado las diferencias dejamos de lado el efecto de mayor
penalización a mayores distancias, y todas las distancias se penalizan
por igual.
De este modo, para minimizar la función, deberíamos
lograr obtener el vector que tenga, en total, las menores distancias.
Para esto necesitamos el valor que se encuentre exáctamente en el
“medio” de los datos, es decir, la mediana.
norma <- function(vector_aleatorio) {
return(sqrt(sum(vector_aleatorio ^ 2)))
}
dame_tita <- function(n_var_al, n_angulos) {
titas = c()
for (i in seq(1, n_angulos)){
X = rnorm(n_var_al)
Y = rnorm(n_var_al)
producto_escalar = X %*% Y
normaX = norma(X)
normaY = norma(Y)
tita_rad = acos(producto_escalar / (normaX * normaY))
titas[i] = tita_rad
}
return (titas);
}
for (n in c(2,3,4,5,10,20,30,100)){
trial = c()
var_trial = c()
for (i in 1:1000) {
trial[i] = mean(dame_tita(n, i))
var_trial[i] = var(dame_tita(n, i))
}
plot(var_trial, col = "green",
main=paste("VARIANZA PARA n = ",n),
type = "l",
xlab = "CANTIDAD DE ANGULOS CALCULADOS",
ylab = "VARIANZA DE TITA")
hist(trial,
main = paste0("DISTRIBUCION DE THETAS PARA n = ", n),
xlab = "VALOR DE THETA [rad]",
ylab = "FRECUENCIA",
col = "blue")
plot(seq(1:1000), trial, type='l',
main=paste('LEY DE GRANDES NUMEROS PARA n = ',n),
xlab = "CANTIDAD DE ANGULOS PROMEDIADOS",
ylab = "ESPERANZA DE THETA",
ylim=c(.5, 3))
abline(h=1.56, col='red')
abline(h=1.56*1.05, col='red', lty=2)
abline(h=1.56*.95, col='red', lty=2)
}
Notamos como la varianza de theta disminuye a medida que aumenta el numero de dimensiones. Esto lo vemos claramente al comparar los graficos de varianza para n = 2 variables aleatorias, donde la misma toma un valor en torno a 0.82, y la varianza para n = 100 variables aleatorias, donde toma un valor en torno a 0.01.
library('plotly')
## Loading required package: ggplot2
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
#Funciones
muestreo_poblacional = function(casos_p, verbose=F){
max_tiempo = 10
tiempo = seq(1:max_tiempo)
#Usamos una matriz para ir guardando los distintos valores que vayan saliendo, al final, cada columna de la matriz será un sample de los 100 que se hacen, es por esto que la matriz será de (10,100) --> (tiempos, iteraciones)
poblacion_final = matrix(ncol=0,nrow=max_tiempo)
for(k in 1:200){
if (verbose==T){
print(paste('Round: ',k))
}
#Comenzamos con 1 solo sobreviviente
poblacion = c(1)
#Vector para ir guardando los valores en cada tiempo t
poblacion_total = c()
for (t in tiempo){
#Se toma un sample de una distribución aleatoria uniforme, si el valor/100 es menor a casos_p, el elemento sobrevive,
#caso contrario el mismo muere
probas = runif(length(poblacion), 1, 100) / 100
poblacion = probas * poblacion
sobrevivientes = poblacion < casos_p
#Lo población es igual a todas las divisiones de los sobrevivientes
poblacion = rep(1, times = sum(sobrevivientes) * 2)
poblacion_total[t] = sum(poblacion)
}
#Populamos la matriz
poblacion_final = cbind(poblacion_final, poblacion_total)
}
return(poblacion_final);
}
crecimiento_conocido = function(t, casos_p){ #Esta es la función conocida, a la que nuestra curva tiene que parecerse
return((2 * casos_p) ^ t)
}
evolucion_poblacional = function(poblacion_final, casos_p){
for (i in c(1,2,3,4,7,20,40,80,95,100,150,200)){
#plot(seq(1,10), apply(poblacion_final, 1, FUN=first), type='l', lwd=2, col='red')
if (i==1){
fig = plot_ly(width = 1500, height = 500, name='Sample') %>%
add_trace(x = seq(1,10), y = poblacion_final[,i], type = 'scatter', mode = 'lines+markers',
marker = list(line = list(width = 2,color='darkgreen')),
name = 'Muestra',
line=list(color='darkgreen')) %>%
add_trace(x=seq(1,10), y=crecimiento_conocido(seq(1,10),casos_p), type = 'scatter', name = 'Teórica', mode = 'lines+markers',
line=list(color='red'), marker=list(color='red'), name='Teórico')
} else {
fig = plot_ly(width = 1500, height = 500, showlegend = FALSE) %>%
add_trace(x = seq(1,10), y = apply(poblacion_final[,1:i], 1, FUN=mean), type = 'scatter', mode = 'lines+markers',
marker = list(line = list(width = 2,color='darkgreen')),
name = 'Muestra',
line=list(color='darkgreen')) %>%
add_trace(x=seq(1,10), y=crecimiento_conocido(seq(1,10), casos_p), type = 'scatter', name = 'Teórica', mode = 'lines+markers',
line=list(color='red'), marker=list(color='red'))
}
assign(
paste0(
"plot_", i),fig)
}
fig <- subplot(plot_1, plot_2, plot_3, plot_4, plot_7, plot_20, plot_40,
plot_80, plot_95, plot_100, plot_150,plot_200,
nrows = 2,
shareY = TRUE,
margin = 0.01) %>%
layout(title = paste('Evolución de la población con LGN vs curva teórica (P=',casos_p,')'),
annotations = list(
list(x = 0.07, y = .9, text = "Iteración 1", xref = "paper", yref = "paper", xanchor = "center", yanchor = "bottom", showarrow = FALSE),
list(x = 0.24, y =.9, text = "Iteración 2", xref = "paper", yref = "paper", xanchor = "center", yanchor = "bottom", showarrow = FALSE),
list(x = .41, y = .9, text = "Iteración 3", xref = "paper", yref = "paper", xanchor = "center", yanchor = "bottom", showarrow = FALSE),
list(x = .58, y = .9, text = "Iteración 4", xref = "paper", yref = "paper", xanchor = "center", yanchor = "bottom", showarrow = FALSE),
list(x = .75, y = .9, text = "Iteración 7", xref = "paper", yref = "paper", xanchor = "center", yanchor = "bottom", showarrow = FALSE),
list(x = .92, y = .9, text = "Iteración 20", xref = "paper", yref = "paper", xanchor = "center", yanchor = "bottom", showarrow = FALSE),
list(x = 0.07, y = .35, text = "Iteración 40", xref = "paper", yref = "paper", xanchor = "center", yanchor = "bottom", showarrow = FALSE),
list(x = 0.24, y =.35, text = "Iteración 80", xref = "paper", yref = "paper", xanchor = "center", yanchor = "bottom", showarrow = FALSE),
list(x = .41, y = .35, text = "Iteración 95", xref = "paper", yref = "paper", xanchor = "center", yanchor = "bottom", showarrow = FALSE),
list(x = .58, y = .35, text = "Iteración 100", xref = "paper", yref = "paper", xanchor = "center", yanchor = "bottom", showarrow = FALSE),
list(x = .75, y = .35, text = "Iteración 150", xref = "paper", yref = "paper", xanchor = "center", yanchor = "bottom", showarrow = FALSE),
list(x = .92, y = .35, text = "Iteración 200", xref = "paper", yref = "paper", xanchor = "center", yanchor = "bottom", showarrow = FALSE)
)
)
return (fig);
}
P=.2
poblacion_final = muestreo_poblacional(.2)
fig = evolucion_poblacional(poblacion_final, .2)
fig
P=.5
poblacion_final = muestreo_poblacional(.5)
fig = evolucion_poblacional(poblacion_final, .5)
fig
p=.7
poblacion_final = muestreo_poblacional(.7)
fig = evolucion_poblacional(poblacion_final, .7)
fig
p=.8
poblacion_final = muestreo_poblacional(.8)
fig = evolucion_poblacional(poblacion_final, .8)
fig
INTERPRETACION: Observamos que a mayor p, la convergencia con la curva teorica se da de manera mas rapida.
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
datos <- readRDS(file="data_final_grupal.Rda")
3a) i) Realice un histograma del peso de los recien nacidos. Hint: hist(⋆)
hist(datos$peso)
hist(datos$peso, xlab = "Peso", main = "Histograma del peso", col=2)
(mean(datos$peso, na.rm=T))
## [1] 3264.672
(median(datos$peso, na.rm = T))
## [1] 3310
sd(datos$peso, na.rm = T)
## [1] 595.2408
mad(datos$peso, na.rm = T)
## [1] 493.7058
IQR(datos$peso, na.rm = T)
## [1] 657
datos$peso %>% boxplot
datos$peso %>% boxplot(col = "blue")
plot(ecdf(datos$peso))
qqnorm(datos$peso)
cuantiles_teoricos = qexp(c(1:length(datos$peso)) / (length(datos$peso) + 1), 1)
qqplot(cuantiles_teoricos, datos$peso,
xlab = "Cuantiles Teoricos Exp(1)",
ylab = "Cuantiles muestrales")
summary(datos$`edad madre`, na.rm = T)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 12.00 23.00 27.00 27.37 32.00 50.00
max(datos$`edad madre`, na.rm = T)
## [1] 50
min(datos$`edad madre`, na.rm = T)
## [1] 12
(mean(datos$`edad madre`, na.rm=T))
## [1] 27.37429
(median(datos$`edad madre`, na.rm = T))
## [1] 27
sd(datos$`edad madre`, na.rm = T)
## [1] 6.164848
mad(datos$`edad madre`, na.rm = T)
## [1] 7.413
IQR(datos$`edad madre`, na.rm = T)
## [1] 9
4.- Variable categórica. (a) Estudie la variable tipo parto.
(a <- table(datos$`tipo parto`))
##
## C-section Unknown Vaginal
## 132068 1517 293738
table(datos$`tipo parto`)[1]
## C-section
## 132068
barplot(a, col = "violet")
pie(a)
datos$dia %>%
table() %>%
barplot(col = "pink")
dias <- c("lunes", "martes", "miercoles", "jueves", "viernes", "sabado","domingo")
datos$dia <- factor(datos$dia, levels = dias)
datos$dia %>%
table() %>%
barplot(col = "pink")
Estudiando la relación entre dos variables
5.- Relación entre una variable numerica y una categorica. (a) Estudie la relacion entre peso del recien nacido y la multiplicidad del parto. i. Grafique esta relacion plot(datos\(parto multiple, datos\)peso)
datos %>%
select(c("parto multiple", "peso")) %>%
plot()
Cuantos más chicos nacen en el parto, menor es el peso de cada bebé por el lugar que ocupan en el vientre, y el tamaño maximo que el mismo se puede expandir.
library("ggplot2")
graf <- datos %>% select(educ, `edad madre`)
ggplot(graf, aes(educ, `edad madre`)) +
geom_boxplot() +
theme(axis.text.x = element_text(angle = 90, hjust=1))
En la mayoria de los casos vemos que las familias se originan entre los 20 y 30 años de la madre. Los mayores niveles de dispersion de edad en madres, los vemos en los casos “null, not on certificate y no formal education” La dispersion de la edad de las madres en los niveles educativos basico y secundario, son similares. La menor dispersion de los años en las madres, las vemos en los casos que estan cursando su 2do y 3er año estudios universitarios, probablemente porque son nacimientos mas planificados, a diferencia de los otros casos que son mas espontaneos. Tambien vemos en estos casos, que la mediana de las madres es superior a los otros escenarios.
elementary <- c("1 Years of elementary school", "2 Years of elementary school",
"3 Years of elementary school", "4 Years of elementary school",
"5 Years of elementary school", "6 Years of elementary school",
"7 Years of elementary school", "8 Years of elementary school")
high_school <- c("1 year of high school", "2 years of high school",
"3 years of high school","4 years of high school")
college <- c("1 year of college", "2 years of college", "3 years of college")
other <- c("Not on certificate", "NULL", "No formal education")
datos <- datos %>%
mutate(nivel_estudio = case_when(educ %in% elementary ~ "elementary",
educ %in% high_school ~ "high school",
educ %in% college ~ "college",
educ %in% other ~ "other"))
graf <- datos %>% select(educ, `edad madre`, nivel_estudio)
ggplot(graf, aes(educ, `edad madre`, color = nivel_estudio)) +
geom_boxplot() +
theme(axis.text.x = element_text(angle = 90, hjust=1))
6.- Relacion entre una variable categorica y otra categorica.
#Estudio de 2 variables categoricas
table(datos$`tipo parto`, datos$dia)
##
## lunes martes miercoles jueves viernes sabado domingo
## C-section 20454 22921 23103 22825 23233 10696 8836
## Unknown 272 247 252 258 289 109 90
## Vaginal 42031 46607 46935 47081 44858 34878 31348
tabla <- table(datos$`tipo parto`,datos$dia)
barplot(tabla, beside = T , col = c("red", "green"))
plot(datos$gestacion, datos$apgar)
No se puede observar ninguna relacion clara
6.- Relacion entre una variable categorica y otra categorica.
aggregate(datos$apgar,
by = list(datos$gestacion),
mean, na.rm=T)
## Group.1 x
## 1 12 9.000000
## 2 15 8.500000
## 3 17 1.529412
## 4 18 2.075000
## 5 19 1.846154
## 6 20 1.823009
## 7 21 1.610390
## 8 22 1.920792
## 9 23 4.115789
## 10 24 5.439276
## 11 25 6.330275
## 12 26 6.664444
## 13 27 7.147482
## 14 28 7.466565
## 15 29 7.667149
## 16 30 7.889899
## 17 31 8.081914
## 18 32 8.260337
## 19 33 8.399011
## 20 34 8.548649
## 21 35 8.662334
## 22 36 8.775058
## 23 37 8.872005
## 24 38 8.926912
## 25 39 8.933183
## 26 40 8.919064
## 27 41 8.891297
## 28 42 8.871996
## 29 43 8.734848
## 30 44 8.906250
## 31 45 8.666667
## 32 46 9.200000
## 33 47 9.000000
## 34 48 9.500000
## 35 51 9.000000
(tabla3 = aggregate(datos$apgar,
by = list(datos$gestacion),
mean, na.rm=T))
## Group.1 x
## 1 12 9.000000
## 2 15 8.500000
## 3 17 1.529412
## 4 18 2.075000
## 5 19 1.846154
## 6 20 1.823009
## 7 21 1.610390
## 8 22 1.920792
## 9 23 4.115789
## 10 24 5.439276
## 11 25 6.330275
## 12 26 6.664444
## 13 27 7.147482
## 14 28 7.466565
## 15 29 7.667149
## 16 30 7.889899
## 17 31 8.081914
## 18 32 8.260337
## 19 33 8.399011
## 20 34 8.548649
## 21 35 8.662334
## 22 36 8.775058
## 23 37 8.872005
## 24 38 8.926912
## 25 39 8.933183
## 26 40 8.919064
## 27 41 8.891297
## 28 42 8.871996
## 29 43 8.734848
## 30 44 8.906250
## 31 45 8.666667
## 32 46 9.200000
## 33 47 9.000000
## 34 48 9.500000
## 35 51 9.000000
barplot(tabla, beside = T , col = c("red", "blue","green"))
legend("topright",
c("C-section", "Unknown", "Vaginal"),
fill = c("red", "blue", "green"))
plot(datos$gestacion, datos$apgar)
Este grafico nos dice muy poco de los datos. Esperaríamos que a mayor cantidad de semana de gestación, mayor sea el índice. Sin embargo, se ven muchos puntos dispersos porque este gráfico no tiene en cuenta la superposición de puntos.
tabla3 = aggregate(datos$apgar,list(datos$gestacion), mean,na.rm=T)
plot(tabla3,xlab="gestacion",ylab="apgar")
Ahora sí encontramos una relación más clara. A medida que aumentan las semanas de gestacion del bebe, aumenta el índice de apgar.